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AN  ALGORITHMIC  ANALYSIS  OF  THE  MMPP/G/ 1  QUEUE 


Levent  Gun* 

Electrical  Engineering  Department  and  Systems  Research  Center 
University  of  Maryland,  College  Park,  Maryland  20742. 


ABSTRACT 

A  single  server  queue  with  general  service  time  distribution  is  considered  when 
the  input  is  a  Markov  modulated  Poisson  process  (MMPP).  An  algorithmic  solution 
to  the  transform  of  the  stationary  delay  and  queue  length  distributions  is  summa¬ 
rized,  and  recursive  closed-form  expressions  are  obtained  for  the  moments  of  these 
distributions.  The  numerical  implementation  of  these  results  is  discussed  in  detail 
with  particular  reference  to  an  algorithm  due  to  Lucantoni  and  Ramaswami  [11]  and 
its  accelerated  version  due  to  Ramaswami  [19].  This  algorithm  is  shown  to  be  an 
efficient  tool  in  the  matrix-analytic  solution  of  many  stochastic  models,  as  various 
steps  for  saving  considerable  amounts  of  unnecessary  computations  are  identified. 
A  special  case  of  the  model  where  the  service  time  distribution  is  of  phase  type  is 
discussed  and  the  stationary  queue  length  distribution  at  arbitrary  times  is  obtained 
in  matrix- geometric  form.  Finally,  the  matrix-geometric  and  the  M /G / 1  approaches 
are  compared  through  this  special  case. 
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1.  INTRODUCTION 

The  basic  motivation  behind  the  work  reported  in  this  paper  is  to  explore 
the  computational  aspects  of  matrix-analytic  solution  techniques  in  the  study  of 
many  probabilistic  models  arising  in  applications.  These  methods  were  pioneered 
by  Neuts  and  have  been  applied  successfully  to  solve  many  stochastic  models  [7,  10- 
15].  Central  to  the  matrix-analytic  approach  is  the  entrywise  minimal  nonnegative 
solution  of  non-linear  matrix  equations  of  the  form 

OO  CO 

G=J2AnGn  or  R  =  Y^RnAn.  (1.1) 

71=0  71=0 

Several  iterative  algorithms  for  solving  these  matrix  equations  are  reviewed  in  the 
Appendix.  For  the  M/G/l  type  queues  such  as  the  PH/G/1  and  the  MMPP/G/ 1 
queues,  an  efficient  and  numerically  stable  algorithm  to  compute  these  matrices 
without  having  to  compute  and  store  the  matrices  {An}§°  was  given  in  the  work  of 
Lucantoni  and  Ramaswami  [11].  In  this  paper,  insights  gained  on  the  behavior  of 
this  algorithm  through  numerical  studies  are  reported.  In  view  of  these  observations, 
a  simple  extension  based  on  linear  extrapolation  is  proposed  for  the  computation  of 
the  G  matrix.  This  minor  extension  leads  to  a  new  stopping  criterion  and  results 
in  considerable  amounts  of  savings  in  the  computations  for  a  given  accuracy. 

This  paper  illustrates  the  above  mentioned  methodology  and  algorithmic  re¬ 
sults  by  studying  both  the  analytical  and  computational  aspects  of  the  stationary 
waiting  time  and  the  queue  length  distributions  for  MMPP/G/ 1  queueing  sys¬ 
tems.  The  MMPP/G /I  queue  is  a  single  server  system  with  general  service  time 
distribution  where  the  arrivals  are  modeled  by  a  Markov  modulated  Poisson  Pro¬ 
cess  (MMPP)',  the  first-come-first-served  (FCFS)  service  discipline  is  enforced.  The 
MMPP  is  a  doubly  stochastic  Poisson  process  whose  rate  is  determined  by  the  state 
of  a  continuous-time  Markov  chain,  and  it  has  great  appeal  in  its  applicability  in 
modeling  many  problems  [6,  7]. 

The  MMPP/G/ 1  queue  is  a  special  case  of  the  more  general  N/G/l  queue 
studied  by  Ramaswami  [17].  Most  of  the  analytical  results  mentioned  here  can 
also  be  obtained  directly  from  the  results  obtained  in  [17]  by  specializing  the  N- 
process  to  a  MMPP.  However,  the  simplicity  of  the  MMPP/G /I  queue  leads  to 


-  1  - 


closed-form  expressions  in  the  recursive  calculation  of  higher  order  moments  of  the 
above-mentioned  distributions. 

This  paper  is  organized  as  follows.  In  Section  2,  the  queueing  model  is  described 
and  equations  describing  the  transition  probability  matrix  are  obtained.  In  Section 
3,  an  algorithm  from  reference  [7]  is  briefly  outlined  to  compute  the  distributions 
of  the  virtual  waiting  time  and  of  the  waiting  time  seen  by  an  arriving  customer. 
In  Section  4,  transforms  of  the  queue  length  distributions  at  departure  epochs  and 
at  arbitrary  times  are  obtained.  In  both  Sections  3  and  4,  recursive  expressions 
are  given  for  all  the  moments  of  these  distributions.  Section  5  contains  various 
insights  obtained  through  numerical  experimentation  of  the  Lucantoni-Ramaswami 
algorithm  on  the  MMPP/G/ 1  queue.  An  accelerated  version  of  this  algorithm  due 
to  Ramaswami  is  also  discussed  and  numerical  comparisons  are  made.  Finally,  in 
Section  6,  the  stationary  queue  length  distribution  and  its  first  two  moments  are 
obtained  at  an  arbitrary  time  in  the  special  case  where  the  service  time  distribution 
is  of  phase  type.  Computational  requirements  in  this  special  case  are  also  discussed 
and  compared  with  the  algorithm  available  for  the  general  service  distributions. 

A  word  on  the  notation  used  hereafter:  The  r  x  r  identity  matrix  is  denoted  by 
Ir  and  the  r  x  1  column  vector  of  ones  is  denoted  by  er,  while  the  1  x  r  dimensional 
row  vector  with  zero  entries  is  denoted  by  0r.  The  notation  O  is  used  for  the 
zero  matrix  with  appropriate  dimensions.  Also,  the  spectral  radius  of  matrix  X  is 
denoted  by  sp(X )  and  “T”  is  used  to  denote  the  transpose  operator. 

2.  THE  MODEL 

Consider  an  m-state  irreducible  continuous-time  Markov  chain  with  infinitesi¬ 
mal  generator  matrix  S  and  the  stationary  probability  distribution  vector  7 r.  When 
the  Markov  chain  is  in  state  i,  1  <  i  <  m,  the  arrivals  are  Poissonian  with  rate 
A i  and  fed  into  a  FCFS  single  server  queue  whose  service  times  are  independent 
and  identically  distributed  with  common  distribution  function  H(-).  Let  A  be  the 
m  x  1  column  vector  with  ith  component  A;.  The  kth  moment  of  the  service  time 
distribution  about  the  origin  is  denoted  by  and  the  m  x  m  diagonal  matrix 
with  elements  A ;  along  the  diagonal  is  denoted  by  A.  It  is  assumed  that  the  system 
parameters  are  such  that  the  queueing  system  is  a  stable  one,  i.  e.,  ir\p^  <  1.  Let 
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{ rn  :  n  >  0}  denote  the  successive  departure  epochs,  with  To  =  0,  and  define  X(t) 
and  J(t )  to  be  the  number  of  customers  in  the  system  and  the  phase  of  the  MMPP 
at  time  t+ ,  respectively.  The  sequence  {(X(r„),  J(r„),  r„+1  —  r„)  :  n  >  0}  form 
a  semi-Markov  sequence  on  the  state  space  {0, 1,. ..}  X  {1, . . . ,m }  with  transition 
probability  matrix  Q(-)  given  by 


/  B0{x) 

Bi(x) 

B2(x)  ■ 

*  \ 

A0(x) 

Ai(x) 

A2(x)  ■ 

O 

A0(x) 

Ai(x)  ■ 

O 

O 

A0(x)  ■ 

V  ! 

‘ 

[  ' 

i  7 

(2.1a) 


where  the  m  X  m  block  entries  are  given  by 


-|4-n(*^)  — 

/  P(n,t)  dH(t)  ,  n  —  0, 1,  •  •  • 

Jo 

,  x  >  0, 

(2.16) 

Bn(x)  = 

U  *  An(x),  n  -  0, 1,  •  -  • 

,  x  >  0, 

(2.1c) 

with 


Pij(n,t )  =  P{X(t)  =  n,  J(t)  =  ,/|X(0)  =  0,  J(0)  =  i}  ,  1  <  i,j  <  rn  ,  (2.1c) 


and 


U(x)=  [  P(0,t)  A 

J  o 


dt . 


(2. Id) 


The  operator  *  in  equation  (2.1c)  is  the  matrix  convolution  operator.  By  using 
standard  techniques  [8],  the  matrices  P(n,t )  can  be  shown  to  satisfy  the  following 
Chapman-Kolmogorov  equations 


dP(n,  t ) 
dt 


P(n ,  t )  (S  —  A)  +  P{n  —  1,  t)  A  ,  if  n  —  1, 2,  • 


P(0,  t)(S~  A), 


if  n  =  0. 


(2.2) 


3.  STATIONARY  WAITING  TIME  DISTRIBUTIONS 

An  algorithm  to  compute  the  distributions  of  the  virtual  waiting  time  and  of  the 
waiting  time  seen  by  an  arriving  customer  is  given  in  [7].  For  sake  of  completeness, 
the  steps  of  this  algorithm  are  briefly  outlined  below. 
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1.  Compute  the  m  x  m  irreducible  stochastic  matrix  G  and  its  stationary  prob¬ 

ability  vector  g.  The  G  matrix  is  the  entry  wise  minimal  nonnegative  solution 
of  the  non-linear  matrix  equation  G  =  AnGn,  with  An  :=  A„(oo).  For 

1  <  hj  <  mi  Gij  denote  the  probability  that  a  busy  period  starting  with  the 
MMPP  in  phase  i  ends  in  phase  j.  The  matrix  G  is  a  key  ingredient  in  ob¬ 
taining  most  quantities  of  interest  such  as  the  busy  period  and  waiting  time 
characteristics  and  is  studied  extensively  by  Neuts  [14].  An  iterative  procedure 
for  computing  the  matrix  G  was  given  by  Lucantoni  and  Ramaswami  [11]  and 
is  discussed  in  Section  5. 

2.  Compute  the  m  x  m  stochastic  matrix  A  :=  —  /o°°  eSt  dH(t).  Here, 

for  1  <  i,j  <  m,  Aij  is  the  probability  that  a  service  time  ends  with  the 
MMPP  in  phase  j  given  that  the  service  began  in  phase  i. 

3.  Compute  the  m  X  m  stochastic  matrix  U  :=  U( oo)  =  (A  —  S’)-1  A.  For  1  < 
i,j  <  m,  Uij  is  the  probability  that  the  first  arrival  to  a  busy  period  arrives 
with  the  MMPP  in  phase  j  given  that  the  last  departure  from  the  previous 
busy  period  departed  with  the  MMPP  in  phase  i. 

4.  Compute  the  m  x  1  column  vectors 

/?  =  /i(1)  (xA)em  +  (S  +  em7r)~1  (A  -  Im) A 

and 

/r  —  (An  G  T  em<7)[Fm  AT  *  @g\  €rn  • 

For  1  <  i  <  m,  the  ith  components  fa  and  \ii  are  the  expected  number  of 
arrivals  during  a  service  that  began  in  phase  i  and  the  expected  number  of 
departures  during  a  busy  period  that  began  in  phase  i,  respectively. 

5.  Compute  the  1  X  m  row  vector  x0  =  (dUg,)~1d,  where  the  mxl  row  vector  d  is 
such  that  dU G  =  d  and  dem  =  1,  with  the  interpretation  that  the  ith  component 
of  xq  is  the  stationary  probability  that  a  departure  leaves  the  system  empty  with 
the  MMPP  in  phase  i. 

6.  Compute  the  lxm  row  vector  yo  =  (x  A)  xo  (A  —  S')-1  where  the  ith  component 
of  2/0  is  the  stationary  probability  of  the  system  being  empty  and  that  the  phase 
of  the  MMPP  is  i  at  an  arbitrary  time. 
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7.  The  Laplace-Stieltjes  transform  (LST)  of  the  virtual  delay  distribution  is  given 
by  W(s)em  where 

_  {  sy0  [slm  +  S  -  A  +  A^(s)]-1  if  5  >  0, 

W{*)  =  {  (3-1) 

f  7T  if  s  =  0, 

and  H(-)  is  the  LST  of  LT(-).  The  ith  component  of  W(x),  the  inverse  transform 
of  the  lxm  row  vector  W{s),  is  the  joint  probability  that  at  an  arbitrary  time 

the  MMPP  is  in  phase  i  and  that  a  virtual  customer  who  arrived  at  that  time 

would  wait  less  than  or  equal  to  x  time  units  before  entering  service. 

8.  Finally,  the  waiting  time  distribution  Wa{-)  seen  as  by  an  arrival  is  given  by 

Wa(x)  =  (7tA)-1  W(x)\  .  (3.2) 


The  first  two  moments  of  the  virtual  waiting  time  distribution  are  available  in 
[7].  Here,  in  view  of  the  equations  (3.1)  and  (3.2)  recursive  expressions  are  given  for 
calculating  arbitrary  moments  of  the  waiting  time  distributions.  Let  p  be  the  traffic 
intensity  given  by  p  —  ttX and  pose 


iy(")  := 


I'F(O) 


£-sw(s) 


s=0+ 


if  n  =  0, 

if  n  =  1,  2,  •  •  •  , 


Lemma  3.1  The  moments  of  the  virtual  waiting  time  distribution  are  given  by 
E{Wn)  =  (-1  )nW(n)em 

x_j\n+i  (3.3  a) 

=  (n  +  1)(1  -  p)^n  +  1^(1)rfn  +  Cn+ ^  A  ’  71  =  1> 2’  ‘  ‘ ' 

where  the  lxm  vectors  {c„}£°  and  are  defined  by 


:=  ^(-l)Vm)(n)w("_m)  , 

m  =  2  \m/ 


dn  . — 


[tt(/to  -  /r(1)A)  -  2/0 ]  (s  +  em7T)  1  , 
n?F()l-1)(Jm  -  /i(1)A)  +  cnA]  (5  +  em7r)-\ 


=  2,3,--- 


if  n  =  1 
if  n  =  2,3,  •  •  •  . 
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The  1  x  m  vectors  {P0n)}g°  in  these  definitions  are  calculated  hy  the  relations 

W{0)  =  7T  ,  W(n)  =  (-lfn)E(Wn)ir-dn  n  =  1,2,  ■  •  •  .  (3.3 b) 

Proof:  Rewrite  equation  (3.1)  as 

W(s)[sl+S-  k  +  KH{s)\  =  sy0  ,  s>0, 

and  differentiate  it  n  times.  Upon  letting  s  — >  0+  in  the  resulting  relation  and  using 
the  definition  of  cn 

_  f  2/o  -  nVm  ~  /r(1)A]  ,  if  n  =  1, 

IU(n)6'=^  _  (3.4) 

[  -nW^~l)  [Im  -  jt<1)  A]  -  cn A  ,  if  n  =  2, 3,  •  •  •  . 

Add  W^emit  to  both  sides  of  equation  (3.4).  Noting  that  n(S  +  emn )  1  =  7r 
and  using  the  definition  of  dn  give  equation  (3.3b)  for  n  =  1,2,  Equation 
(3.3b)  is  trivially  true  for  n  =  0  since  bU(0)  =  n.  Now  replace  n  with  n  +  1  in 
(3.4)  and  premultiply  by  em,  the  desired  result  (3.3a)  is  readily  obtained  after  some 
simplifications.  □ 

The  moments  of  the  stationary  waiting  time  distribution  as  seen  by  an  arbitrary 
arrival  now  follows  from  equation  (3.2)  and  are  given  by 

E(WZ)  =  (-l)n(Tr\)-1W^  \  ,  n  =  1,2,....  (3.5) 

4.  THE  STATIONARY  QUEUE  LENGTH  DISTRIBUTIONS 

In  this  section,  results  from  the  reference  [17]  are  specialized  to  obtain  equations 
satisfied  by  the  transforms  of  the  stationary  queue  length  distributions  at  points  of 
departures  and  at  arbitrary  times.  These  equations  are  then  used  to  generate  closed 
form  recursive  expressions  to  obtain  the  moments  of  these  distributions. 

Let  the  stationary  queue  length  density  at  a  departure  point  be  denoted  by  the 
row  vector  x  which  is  partitioned  into  1  X  m  row  vectors  as  x  :=  (xo,  x\, . . .),  where 
the  ith,  1  <  i  <  m,  component  of  x k  is  the  joint  probability  that  at  a  departure 
epoch  there  are  k  customers  in  the  system  and  that  the  MMPP  is  in  phase  i.  The 
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vector  x  is  the  invariant  probability  vector  of  the  irreducible  stochastic  matrix  Q( oo) 
and  thus  satisfies  the  equations 

x  Q(oo)  =  x  ,  2:600  =  1.  (4.1) 

Equation  (4.1)  can  be  rewritten  in  terms  of  the  block  entries  as 

*+1 

x k  —  £0  Bk  4“  'y  %i  Ak— ,  k  —  0, 1,  •  •  •  ,  (4.2) 

1=1 

where  Bk  :=  Bk( 00),  k  =  0, 1,  •  •  •.  Multiplication  of  equation  (4.2)  by  the  complex 
number  zk  and  summing  over  A;  =  0,1,---,  yield  the  equation 

X(z)  [zlm  -  A{z)\  =  x0  [zU  -  Im\  A(z)  ,  \z\  <  1,  (4.3) 

where 

OO  OO 

X(z)  ^2  %k  zk  and  A(z)  :=  ^  zk  ,  |z|  <  1 . 

fc=0  k=0 

From  equation  (2.1b)  it  easily  follows  that 

r  OO 

A(z)  =  /  P(z,t)dH(t ),  (4.4) 

Jo 

where 

OO 

P(z,t):=  J2  P{k,t)zk  ,  |z|  <1. 

k=0. 

In  view  of  equation  (2.2),  P(z,t )  satisfies  the  following  differential  equation 

^  —  =  P(z,  t )  (S  -  A)  +  zP(z,  t) A  ,  t  >  0 

P(z,  0)  =  Im  , 

which  is  readily  solved  to  yield 

P(z,t)  =  5  |^|  <  1  and  t  >  0.  (4. 5) 
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Therefore,  the  matrix  A(z)  of  equation  (4.3)  is  given  by 


/*oo 

A(z)  =  /  dH(^  ?  \z\<l  . 

Jo 


(4.6) 


Next,  an  explicit  expression  is  obtained  for  the  quantity  X(l~)  =  Xk' 

which  is  needed  in  the  computation  of  the  moments  of  the  queue  length  distribution. 
Letting  1—  in  equation  (4.3)  yields 

X(l-)  [Im  -  A]  —  x0[U  -  Im]  A  .  (4.7) 


Note  by  the  definition  of  the  stochastic  matrix  A  and  the  row  vector  n  that  7 r  is 
also  the  stationary  probability  vector  of  A.  Addition  of  n  =  (X(l~)em)  7 r  to  both 
sides  of  the  equation  (4.7)  readily  yields 

X(l-)  =  x0(U  -  Im)A(Im  -  A  +  em tt)_1  +  7r  , 


through  the  use  of  the  fact  that  7r(Jm  —  A  +  em7r)_1  =  n.  Invertibility  of  the 
stochastic  matrix  (Im  —  A  +  em7r)  is  shown  in  Kemeny  and  Snell,  [9,  p.  100]. 

Simplicity  of  the  equation  (4.3)  allows  a  recursive  computation  of  the  moments 
of  the  queue  length  distribution  at  departure  epochs.  To  that  end,  pose 


XW  :=X(1"), 


:=  A, 


72  —  0 


x(n)  ;=  J^A\z=1_ ,  :=  ,  n  =  1, 2,  •  •  •  . 

The  results  from  the  next  lemma  can  be  proved  by  using  arguments  similar  to  the 
ones  in  the  proof  of  Lemma  3.1. 


Lemma  4.1  The  moments  of  the  queue  length  distribution  at  departure  epochs  are 
given  by 


E(Xn)  =  X(n)< 


(4.8a) 


(n  +  l)(l-p) 


[(n  +  1  )dnA(1)  +  x0Vn+i  -  cn+i]e 


n  —  1,  2,  ■ 


where  the  m  X  m  matrices  and  the  1  x  m  vectors  and  are 
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defined  by 


Vn  r- 


C-n  • —  \ 


dn  .  - 


(U  -7m)A(n)  +  nUA(n~V  ,  n  =  1, 2,  •  •  •  , 

(  0  if  »  =  1, 

i  E”=2  C)^(B-m)^(Bl)  if  n  =  2, 3,  •  •  •  , 

[*oK  +  cn  -  «A'<n"1>(/m  -  ^(1))](/m  -  A  +  e^)-1 ,  n  =  1,2,  • 


The  1  X  m  vectors  {X^jg0  in  these  definitions  are  calculated  as 


X(0):=X(1“),  X(n>  =  E(Xn)ir  +  dn  ,  »  =  1,2,---.  (4.86) 


Let  the  stationary  queue  length  density  at  an  arbitrary  time  be  denoted  by  the 
row  vector  y  which  is  also  partitioned  into  1  X  m  row  vectors  as  y  :=  (t/g,  yi, . . .), 
where 


(; yk)i~  lim  P[X(t)  =  k,  J (t)  =  i\X (0)  =  k' ,  J(0)  =  i'],  1  <  i,i'  <  m  , 
<—►00 


k,  k'  =  0, 1,  •  •  •  . 


The  following  result  is  a  special  case  of  Theorem  3.3.18  of  [17,  p.  246]  and  can  be 
obtained  using  the  Key  Renewal  Theorem  [2]. 

Lemma  4.2  The  generating  function  Y(z )  =  is  given  by 

(  (7tA)_1(1  -  z)X(z)  [(1  —  z) A  —  5]_1  ,  if  0  <  2  <  1 , 

Y(z)  =  \  (4.9) 

l  7T  ,  if  2  =  1  . 


Use  of  equation  (4.9)  and  calculations  similar  to  the  ones  made  in  the  proof  of 
Lemma  3.1  give  the  following  result. 


Lemma  4.3  The  moments  of  the  stationary  queue  length  distribution  at  arbitrary 
times  are  given  by 


E(Yn)  =  y(")em  ,  n  =  1,2,  •  •  ■  ,  (4.10a) 

=  E(Xn )  -  n  -  (7rA)-1K(n“1)A]  (5  +  em7r)_1A  , 
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where  the  1  x  m  vectors  :=  d  -Jn^  \z=i- ,  ^  =  0,1,---,  are  given  by 


{7T  if  n  =  0, 

n  [(TrA)^^-1)  -  y(ri-1)A]  (5  +  e^)-1  +  E(Yn)jr,  if  n  =  1, 2,  •  •  -. 

(4.106) 

5.  IMPLEMENTATION  DETAILS  -  EXTENSIONS  AND  COMPAR¬ 
ISONS 

The  central  item  in  the  algorithm  outlined  in  Section  3  is  the  so-called  G  matrix 
introduced  by  Neuts  [14].  An  efficient  and  numerically  stable  way  of  computing  the 
G  matrix,  without  having  to  evaluate  the  integrals  in  (3.1),  is  given  in  the  work 
of  Lucantoni  and  Ramaswami  [11].  This  algorithm  is  based  on  the  randomization 
technique  of  Grassman  [3]  and  can  be  applied  to  the  matrix-analytic  solution  of 
many  stochastic  models.  A  sub-Newton  scheme  is  later  incorporated  to  this  algo¬ 
rithm  by  Ramaswami  [19]  to  improve  its  rate  of  convergence.  In  this  section  an 
extrapolation  method  is  suggested  for  these  algorithms  based  on  the  insights  gained 
through  numerical  studies.  Although  this  proposed  extension  is  a  minor  addition  to 
the  existing  algorithms,  it  has  been  observed  to  improve  their  efficiency  considerably 
especially  when  their  convergence  rates  are  slow. 

THE  RANDOMIZATION  ALGORITHM 

This  above  mentioned  algorithm  of  Lucantoni  and  Ramaswami  is  called  the 
Randomization  Algorithm  (RA)  in  the  rest  of  this  paper  and  uses  the  following 
theorem  in  Qinlar  [2]. 

Theorem  5.1  Let  S  be  the  infinitesimal  generator  of  an  m-state  Markov  process 
and  suppose  that  —Su<0<  oc  for  1  <  i  <  m.  Then  the  transition  function 
P(t )  =  exp(St)  may  be  written  as 


n= 0 


where  K  is  the  stochastic  matrix  Im  +  0  1S. 

For  the  model  described  in  Section  2,  the  RA  can  be  given  as 


(5.1) 
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H[n+1)  =  [Im  +  6~1(S -A)]H(kn) +  9-1AHikn)Gk,  n  =  0,1,...,  (5.2) 


Gk+ 1  =  ’  ^  =  0,1,-.., 

n— 0 

where  Co  =  O,  II ^  =  Jm,  A;  =  0, 1,  -  -  - ,  6  —  maxi<j<m(A;  -  So)  and 
7„  =  J0°°  e~6x~£  dll(x),  n  —  0,1,  •  •  ■ .  It  is  shown  in  [11]  that  the  sequence  Gk 
converges  monotonically  to  G.  Note  that  is  the  probability  that  a  service  time 
has  n  epochs  of  a  Poisson  process  with  rate  6.  In  some  useful  special  cases  where 
the  service  distribution  is  either  phase  type  or  deterministic,  7„  may  be  computed 
recursively  without  numerical  integration.  These  two  cases  are  mentioned  below. 

1.  When  H(‘)  is  deterministic  with  mass  a,  then 

7o  =  e~8a  ,  7n  =  — 7n— l  n  =  1,2,  •  •  ■  .  (5.3) 

n 

2.  When  H(-)  is  of  phase  type  with  representation  (a,  A),  it  is  shown  by  Neuts  [14, 
pp.  59,  60]  that  {7n}o°  has  discrete  phase  density  with  representation  (/3,B) 
given  by 

0  =  ea(6Im  -  A)-1  ,  B  =  0(0 Im  -  A)~l  (5.4a) 

and 


fim+i  =  <Xm+\  +  <*(0Im  -  A)-1  A0  ,  B°  =  (0Im  -  A)-1  A0  .  (5.46) 

For  the  notation  and  the  proof  of  this  result,  the  reader  is  referred  to  [14].  The 
probabilities  {7^}^°  can  then  be  computed  by  the  recursion 

70  =  am+i  +  rj(l)A°  ,  7„  =  r}(n+  1)A°,  n=  1, 2,  -  *  *  .  (5.5a) 

where 

r](0)  =  6~1a,  Tj(n  +  1)  =  77(a)  B  ,  n  =  0, 1,  -  -  -  .  (5.56) 

In  practice  the  index  n  in  iteration  (5.2)  is  truncated  at  some  positive  integer 
value  N.  In  view  of  the  probabilistic  interpretation  of  7„  the  truncation  index  N 
can  be  chosen  to  satisfy  7«  <  C- 
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MODIFIED  NEWTON-KANTOROVICH  SCHEME  (MNK) 

The  following  modification  is  proposed  by  Ramaswami  [19]  to  improve  the  rate 
of  convergence  of  the  RA.  For  this  scheme,  called  RA_MNK  hereafter,  the  two 
auxiliary  substochastic  matrices  A;,  i  =  1,2,  are  needed.  These  matrices  are  given 
in  [11]  by 

OO 

Ai  =  nK\U)  , 

n—i 

where  the  matrices  I(\n^  are  recursively  defined  by  Kq0>  =  /,  =0,  and 

i4n+1)  =  x^li  +  e-'is-  A)], 

(5.6) 

K(n+1)  _  R(n)y  +  0-1(5  _  A)]  +  J{W  0-1 A  ?  n  >  o,  i  =  1,2. 

Having  computed  Gk+\  from  the  RA  given  in  equation  (5.2),  the  MNK  step  is  now 
obtained  by  setting 

Zk+1  =  (I-A1)~\Gk+1-Gk)  (5.7  a) 

and  updating  Gk+ 1  by 

Gk+ 1  +  AiZk+i  +  A2(GkZk+i  +  Zk+\Gk )  — >  GT+i  .  (5.7b) 

The  RA_MNK  scheme  then  proceeds  by  using  this  new  Gk+ 1  in  equation  (5.2)  to 
get  the  next  iterate  Gk+  2-  It  is  shown  in  [19]  that  this  scheme  also  converges  to  the 
G  matrix. 

PROPOSED  EXTENSION  BASED  ON  LINEAR  EXTRAPOLATIONS 

The  algorithms  RA  and  RA.MNK  are  stopped  when  all  the  entries  in  the 
successive  iterates  of  the  G  matrix  are  within  a  small  number,  say  e,  of  each  other, 
i.  e., 

max  | (Gk+1)ij  -  (Gk)ij\  <  e  .  (5.8) 

1  <z,j  <.m 

Computational  experience  has  shown  that  when  the  convergence  rate  is  slow  (when 
p  is  large),  even  when  the  successive  iterates  satisfy  (5.8),  say  at  iteration  K,  the 
matrix  Gk  can  be  far  away  from  the  limiting  G  matrix,  i.  e.,  Gij  —  (GK)ij  can  be 
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as  much  as  102  to  104  times  e  for  some  i  and  j  depending  on  p.  However,  the  rows 
of  Gk,  when  considered  as  vectors,  were  observed  to  be  almost  colinear  with  the 
corresponding  rows  of  the  G  matrix.  Therefore,  the  most  recent  iterates,  Gk  and 
Gk- i,  can  be  used  to  obtain  an  approximation  G*K  to  the  stochastic  matrix  G  by 
linear  extrapolation,  i.  e., 


G*k  =  Gk  +  L  [Gk  -  GK- 1]  ,  (5.9) 

where  L  is  a  diagonal  matrix  so  that  the  matrix  G*K  is  stochastic.  Computational 
experience  has  shown  that  even  when  Gk  is  far  away  (relative  to  c)  from  being  a 
stochastic  matrix,  G*K  is  typically  much  closer  to  G  and  use  of  G*K  leads  to  much 
more  accurate  results  in  the  calculation  of  the  performance  measures  of  interest. 

The  above  observations  suggest  that  if  G k  is  the  stochastic  matrix  obtained 
from  the  linear  extrapolation  of  Gk  and  Gk- 1,  then  for  the  same  e,  the  stopping 
criterion 

max  \{Gl+l)lj-(Gl)ii\<t,  (5.10) 

should  be  satisfied  at  an  iteration  K*  which  is  much  smaller  than  K,  and  G*K, 
should  be  closer  to  G  than  Gk-  Indeed,  the  stochastic  matrix  G*K.  were  observed 
to  be  much  closer  to  G  than  Gk  in  all  the  numerical  examples  performed,  and  the 
number  of  iterations  K*  was  much  smaller  than  K  especially  for  large  values  of  p. 
Note  that  the  computation  of  the  matrices  G*k  from  Gk  and  Gk- 1,  k  >  1,  requires 
only  of  order  to2  additional  operations,  whence  do  not  have  a  significant  effect  on 
the  total  CPU  time,  especially  for  large  to,  since  each  iteration  is  of  order  N to3. 
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NUMERICAL  EXAMPLES 


These  claims  are  supported  through  many  numerical  examples,  some  of  which 
are  included  in  Table  5.1,  where  the  RA  and  the  RA_MNK  schemes  are  also  com¬ 
pared  with  and  without  the  proposed  extrapolations.  The  numerical  examples  in 
Table  5.1  correspond  to  the  MMPP/D/1  queue  where  the  service  times  are  chosen 
as  4msec,  6msec,  8msec,  10msec,  ll|msec  and  3jmsec  for  Examples  5. 1.1-5. 1.6, 
respectively.  The  input  process  is  the  superposition  of  m  —  1  identical  two-state 
Markov  processes  with  states  denoted  by  1  and  0  and  the  transition  rate  from  state 
i  to  state  j  is  denoted  by  i,j  =  0, 1.  When  the  process  is  in  state  i  it  produces 
Poisson  arrivals  with  rate  A%  i  —  0, 1.  The  input  process  is  therefore  an  m  state 
MMPP  with  a  tridiagonal  generator  matrix  S  and  arrival  rates  A ,•  given  by 

Si,i— 1  =  (1  l)riO  i  Si,i  =  —  S , 

=  {jn  —  *)r0i  ,  A ,  =  A1(i  —  1)  +  A °(m  —  i) ,  for  1  <  i  <  m . 
Note  that  the  state  of  the  MMPP  denotes  the  number  of  the  component  processes 
that  are  in  state  1,  with  states  1  and  m  of  the  MMPP  corresponding  to  the  situation 
where  all  the  Markov  processes  are  in  state  0  and  in  state  1,  respectively.  The 
following  numerical  values  are  used  in  all  there  examples:  r0i  =  2.5618,  rio  = 
1.586,  A0  =  0.00349,  A1  =  57.416. 

In  Table  5.1,  K  and  K*  denote  the  number  of  iterations  required  to  satisfy 
the  stopping  criterions  (5.8)  and  (5.10),  respectively.  The  expected  value  of  the 
virtual  waiting  time  is  displayed  for  comparison.  The  other  computed  performance 
measures  have  shown  a  typical  behavior.  Here,  EW  and  EW*  correspond  to  the 
expected  virtual  waiting  time  obtained  using  the  stopping  criterions  (5.8)  and  (5.10), 
respectively. 


-14- 


TABLE  5.1 


Example  5.1.1. 

m 

5. 

P  = 

0.351302, 

e7  =  10 

7  (N  = 

9). 

> 

I< 

EW 

CPU 

I\* 

EW* 

CPU 

RA 

i(r3 

20 

3.73625e-3 

1.9 

14 

1.68891e-3 

1.4 

io-5 

42 

1.67364e-3 

3.5 

30 

1.65542e-3 

2.5 

io-7 

64 

1.65464e-3 

4.9 

50 

1.65447e-3 

3.8 

io-9 

85 

1.65446e-3 

6.3 

70 

1.65445e-3 

5.1 

RAJVINK 

10“3 

15 

3.00642e-3 

1.9 

11 

1.67773e-3 

1.6 

io-5 

29 

1.66853e-3 

3.0 

22 

1.65499e-3 

2.4 

io-7 

44 

1.65457e-3 

4.2 

35 

1.65446e-3 

3.6 

io-9 

58 

1.65446e-3 

5.5 

48 

1.65445e-3 

4.4 

Example  5.1.2. 

m  = 

5. 

P  = 

0.526953, 

1 

o 

II 

II 

t- 

11). 

e 

K 

EW 

CPU 

I<* 

EW* 

CPU 

RA 

<n 

1 

o 

49 

1.72419e-2 

4.6 

29 

1.26234e-2 

2.8 

10~5 

109 

1.25098e-2 

9.4 

70 

1.24573e-2 

6.1 

10-7 

171 

1.24560e-2 

14.3 

112 

1.24557e-2 

9.2 

RAJVINK 

10“3 

36 

1.56773e-2 

4.1 

22 

1.25645e-2 

2.7 

10-5 

76 

1.24871e-2 

7.7 

49 

1.24567e-2 

5.0 

10“7 

116 

1.24560e-2 

11.5 

77 

1.24557e-2 

7.6 

Example  5.1.3.  m  —  5.  p  =  0.7026045,  £7  =  10  7  (A  =  13). 


€ 

I< 

EW 

CPU 

K* 

EW* 

CPU 

RA 

io-3 

78 

9.56070e-2 

8.2 

38 

8.10129e-2 

4.2 

io-4 

138 

8.22618e-2 

14.0 

60 

8.07974e-2 

6.1 

io-5 

198 

8.09259e-2 

20.2 

82 

8.07759e-2 

8.3 

IO”7 

320 

8.07750e-2 

32.2 

126 

8.07735e-2 

12.4 

RA_  MNK 

io~3 

59 

9.04182e-2 

6.7 

29 

8.09121e-2 

3.7 

io-4 

99 

8.17241e-2 

10.8 

43 

8.07884e-2 

4.9 

IO-5 

139 

8.08685e-2 

14.8 

58 

8.07749e-2 

6.5 

10-7 

219 

8.07745e-2 

23.0 

86 

8.07735e-2 

9.2 

Example  5.1.4. 

m  =  5. 

p  =  0.8782556,  e7  =  10 

"7  (N 

=  14). 

e 

K 

EW 

CPU 

I\* 

EW* 

CPU 

RA 

IO-3 

115 

0.558680 

11.5 

35 

0.515955 

4.1 

10~4 

261 

0.520344 

24.6 

54 

0.515631 

6.0 

IO"5 

422 

0.516074 

45.1 

72 

0.515603 

7.7 

RAJVINK 

IO"3 

93 

0.544938 

12.1 

27 

0.515794 

3.9 

10“4 

194 

0.518730 

24.3 

39 

0.515619 

5.1 

10~5 

302 

0.515912 

37.0 

51 

0.515601 

6.8 

Example  5.1.5. 

m  =  5. 

p  =  0.995356, 

e7  =  10_ 

7  (N  = 

=  15). 

RA 

10~3 

127 

21.9231 

14.3 

32 

21.8105 

3.9 

io~4 

475 

21.8414 

52.5 

47 

21.8102 

5.5 

io-5 

1539 

21.8175 

170.4 

62 

21.8102 

7.1 

RAJVINK 

IO"3 

111 

21.8998 

15.0 

24 

21.8104 

3.5 

IO"4 

396 

21.8349 

50.5 

35 

21.8102 

5.0 

io-5 

1255 

21.8157 

155.8 

45 

21.8102 

6.2 

Example  5.1.6. 

m  =  12, 

p  =  0.805067 

,  e7  =  10 

-7  (N 

=  13). 

RA 

IO"3 

132 

7.86920e-2 

118.4 

49 

5.91192e-2 

45.6 

io-4 

250 

6.03862e-2 

220.0 

93 

5.85371e-2 

84.9 

10“5 

367 

5.86599e-2 

327.8 

138 

5.84762e-2 

120.5 

RAJVINK 

CO 

1 

o 

102 

7.06491e-2 

106.9 

39 

5.88454e-2 

45.6 

10-4 

178 

5.96593e-2 

187.0 

68 

5.85077e-2 

74.1 

10-5 

254 

5.85873e-2 

264.9 

97 

5.84731e-2 

102.2 
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The  following  observations  immediately  follow; 

•  For  a  given  e,  K*  is  about  35-95%  less  than  K  in  Table  5.1  both  for  the  RA 
and  the  RA-MNK  schemes.  Note  that  the  extrapolations  resulted  in  savings 
both  in  the  number  of  iterations  and  in  the  CPU  times  for  all  values  of  p.  The 
amount  of  savings  increased  with  p;  For  the  low  utilization  of  Example  5.1.1 
savings  of  approximately  35-40%  were  obtained,  while  savings  by  a  factor  of  up 
to  25  were  possible  for  the  very  high  utilization  of  Example  5.1.5,  both  in  RA 
and  RA-MNK. 

More  importantly,  the  performance  measures  obtained  either  from  RA  or 
RA-MNK  without  using  the  extrapolations  were  inaccurate  for  larger  values 
of  e  (10~3,  10~4  and  even  10-5).  On  the  other  hand,  the  extrapolations  re¬ 
sulted  in  much  more  accurate  results  in  the  performance  measures  for  the  same 
e  values.  It  is  clearly  seen  from  Table  5.1  that,  generally  EW*  were  accurate 
up  to  n  digits  after  the  decimal,  where  n  =  log{e).  The  same  behavior  can 
also  be  seen  for  the  average  queue  lengths  in  the  last  column  of  Table  6.1.  The 
value  of  e  used  with  the  new  stopping  criterion  thus  provide  direct  information 
on  the  accuracy  obtained  in  the  calculation  of  the  performance  measures.  On 
the  other  hand,  the  stopping  criterion  (5.8)  gave  no  direct  information  on  the 
accuracy  desired  in  the  calculation  of  the  performance  measures.  Therefore, 
the  use  of  the  stopping  criterion  (5.10)  also  saves  considerable  amounts  of  CPU 
time  by  avoiding  unnecessary  iterations  in  order  to  obtain  a  desired  accuracy. 

•  The  algorithm  RA_MNK  required  approximately  15-25%  fewer  number  of  it¬ 
erations  and  0-15%  less  CPU  time  than  the  RA,  both  with  and  without  the 
extrapolations.  Note  that  although  the  RA_MNK  scheme  was  always  faster 
than  the  RA,  in  some  cases  the  corresponding  CPU  times  were  comparable  due 
to  the  additional  computations  in  the  MNK  step.  The  RA_MNK  scheme  also 
resulted  in  slightly  more  accurate  performance  measures. 

FURTHER  COMMENTS 

Next,  the  effect  of  the  choice  of  e7  is  considered.  In  Table  5.2  ea  is  set  to  10“5 
and  EW*  is  displayed  for  different  values  of  e7  for  the  Example  5.1.4  (resp.  Example 
5.1.5).  The  entries  in  the  table  are  obtained  from  the  RA  (used  with  the  proposed 
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extrapolations),  however  similar  results  are  also  observed  for  the  RA_MNK.  The 
probabilities  7„,  1  <  n  <  N,  are  normalized  to  sum  to  one.  This  normalization 
seems  to  improve  the  accuracy  in  the  obtained  performance  measures  especially  for 
large  e7. 


TABLE  5.2 


Example  5.2.1 

:  m 

=  5, 

p  =  0.526953 

(g  —  10 

e7 

K* 

N 

EW* 

10“2 

69 

5 

0.013167 

10~3 

70 

6 

0.012622 

10-4 

70 

8 

0.012463 

10“5 

70 

9 

0.012458 

10~6 

70 

10 

0.012457 

10"7 

70 

11 

0.012457 

Example  5.2.2 

:  m 

=  5, 

p  =  0.8782556 

=  10 

e7 

I{* 

N 

EW* 

i — i 

0 

1 

to 

73 

7 

0.516296 

10~3 

72 

8 

0.515802 

10-4 

72 

10 

0.515614 

10~5 

72 

11 

0.515605 

10~6 

72 

13 

0.515603 

10-7 

72 

14 

0.515603 

The  figures  in  Table  5.2  indicate  that,  although  the  performance  slightly  de¬ 
grades  as  e7  is  increased,  the  “inaccuracy”  in  calculating  the  (7n}o°  seems  to  have  a 
surprisingly  small  effect  and  up  to  50%  savings  can  be  obtained  in  the  computation 
time  in  cases  when  precision  is  not  of  primary  concern.  Also,  comparison  of  the  last 
two  rows  in  both  examples  shows  that  savings  of  20-25%  can  be  obtained  without 
compromising  on  the  accuracy. 

Considerable  savings  in  computation  time  can  also  be  gained  in  cases  where  the 
S  matrix  is  sparse.  However,  the  computation  time  per  iteration  will  still  approxi¬ 
mately  grow  with  m3  since  the  matrices  H and  Gk  in  equation  (5.3)  will  not  be 
sparse  even  when  the  matrix  S  is  sparse. 

Finally,  the  computation  of  the  matrices  A ^  :=  d  U=i-  ,  k  =  0, 1,  •  •  •,  used 
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in  equation  (4.8),  is  discussed.  Equation  (4.6)  and  Theorem  5.1  yield 


r  oo 

A(k)  =  /  tk  est  dH(t)  hk  , 
Jo 

oo 

=  £7 (nk)Kn  Afc  , 

n=0 


fc  =  0,1, 


(5.11) 


with 

tk  e~et  dH(t)  ,  n  =  0,1,---,  and  k  =  0,  l,**-,  (5.12) 

n! 

where  9  and  the  matrix  K  are  given  in  Theorem  5.1.  It  is  plain  that 

OO 

=  *  =  0,1,  -  .  (5.13) 

n— 0 

In  view  of  equation  (5.11),  once  =  7!°^,  n  =  0, 1, •  •  • ,  is  calculated  either  by 
numerical  integration  or  through  the  recursions  given  by  equations  (5.4)  and  (5.5), 
7 nk\  n  =  0, 1,  •  •  • ,  A:  =  1,2,***,  can  be  obtained  recursively  by 

ln]  =  yy  »  n  =  0, 1,  •  •  •  and  k  =  1, 2,  •  •  •  ,  (5.14) 

and  equation  (5.13)  can  be  used  as  a  truncation  criteria. 

6.  SPECIAL  CASE:  THE  MMPP/PH/l  QUEUE 

In  this  section  the  service  time  distribution  is  specialized  to  be  of  phase  type. 
The  stationary  queue  length  distribution  at  arbitrary  times  is  obtained  in  matrix- 
geometric  form,  and  simple  expressions  are  given  for  the  first  two  moments  of  this 
distribution.  The  stationary  waiting  time  distribution  for  GI/PH/1  queue  was 
derived  in  the  work  of  Ramaswami  and  Lucantoni,  [18,  Thm.  1],  and  will  not  be 
repeated  here. 

Let  the  service  distribution  have  irreducible  phase  representation  (o,  Q )  of  order 
l,  where  a  is  the  lx/  row  vector  of  initialization  probabilities  and  Q  is  the  /  X  l 
matrix  of  transition  rates  among  the  transient  phases  of  the  phase  distribution. 
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Also,  denote  the  I  X  1  column  vector  of  absorption  rates  to  the  absorbing  phase  by 

Q°- 

A  natural  state  space  E  for  this  system  is  given  by 


E  = 


{ (o,0 


{  ( k,i,j ) 


if  k  =  0,  1  <  i  <  m, 

if  k  =  1, 2,  •  •  •  ,  1  <  i  <  m,  1  <  j  <  l. 


where  k  indicates  the  buffer  size,  while  i  and  j  represent  the  state  of  the  M M PP 
and  the  phase  of  the  server,  respectively.  The  service  phase  is  not  defined  when 
the  queue  is  empty.  Since  two  events  cannot  occur  simultaneously  (with  positive 
probability),  it  is  an  easy  exercise  to  show  that  the  underlying  continuous-time 
Markov  chain  is  irreducible. 

The  stationary  queue  length  density  at  an  arbitrary  time  is  again  denoted  by 
the  row  vector  y  which  is  partitioned  as  y  :=  (yo,yi, . . .),  where  y0  is  1  X  m  and  yk, 
k  =  1, 2,  •  •  • ,  are  1  xlm  vectors  and  correspond  to  the  stationary  probabilities  of  the 
states  in  lexicographical  order.  The  corresponding  infinitesimal  generator  matrix  P 
then  takes  the  form 


(  S  —  A  A®  a  0  0 

Im®Q°  Ax  A0  0 


0  0  -\ 

0  0  • 


0  Ai  A\  Aq  0  0 

0  0  A2  A\  Ao  0 


where  the  Im  x  Im  matrices  Ak,  0  <  k  <  2,  are  given  by 

A0  =  A®Ih  Ax  =  (S-A)®Q,  A2  =  Im®Q0a. 


(6.1) 


Here  ®  and  ©  denote  the  Kronecker  product  and  the  Kronecker  sum,  respectively, 

[1]. 

In  order  to  obtain  the  stationary  queue  length  distribution  in  matrix-geometric 
form,  let  the  Im  X  Im  matrix  R  be  the  minimal  nonnegative  solution  to  the  matrix- 
quadratic  equation 

R“  A2  +  R  A\  +  A0  =  O  ,  (6-2) 
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First,  for  a  stable  system,  all  the  eigenvalues  of  R  are  shown  to  be  inside  the  open 
unit  disc. 

Lemma  6.1  The  MMPP/PH/1  queueing  system  is  stable  if  and  only  if  sp(R )  <  1. 

Proof:  Since  P  :=  Ao  +  A\  +  A2  =  5©  (Q  +  Q°a ),  the  matrix  P  is  irreducible  owing 
to  the  irreducibility  of  the  matrices  S  and  Q  -f  Q°a.  Therefore  sp(R)  <  1  if  and 
only  if  pAoeim  <  pA2eim,  where  the  1  X  Im  row  vector  p  is  the  stationary  vector 
of  P  [14,  p.  83].  From  the  structure  of  P  it  is  easy  to  see  that  p  =  tt  0  q  where 
the  1  X  I  row  vector  q  is  the  stationary  probability  vector  of  the  matrix  ( Q  -f  Q°a). 
Therefore,  the  statement  of  the  Lemma  follows  since 

pA2e[m  =  qQ°  —  effective  service  rate 


and  similarly 


p  Aq  eim  =  7r  A  =  e ffective  arrival  rate  . 


n 


/  S- A 

Lemma  6.2  The  matrix  Z  := 

\Im®Q° 


A  CD  a 
A\  T  RA2 


is  an  irreducible  generator 


matrix. 

Proof:  Post  multiplication  of  Z  by  after  some  simplifications  yield 


^  lm+m )  — 


01 

1  Q°) 


ei 


But  since 


®Tm  —  (R2 A2  +  RAi  +  Ao)eim  =  (J;TO  —  R)[ A  0  ej  —  R(em  0  Q°)\ 

and  sp(R )  <  1,  the  equation  [A  0  e;  —  R(em  ®  Q0)]  =  0 Jm  holds.  Therefore,  each 
row  of  the  matrix  Z  sum  to  zero.  On  the  other  hand,  by  the  definitions  of  the  block 
entries  of  the  matrix  Z,  its  off-diagonal  entries  are  all  non-negative  and  the  matrix  Z 
is  a  generator  matrix.  Irreducibility  of  the  Z  matrix  follows  from  the  irreducibility 
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of  the  phase  representation  (a,  Q )  and  of  the  matrix  S,  and  from  the  non-negativity 
of  the  matrix  R.  n 

The  following  matrix- geometric  form  of  the  stationary  probability  vector  y 
now  follows  from  Theorem  1.5.1.  of  [14,  p.  25]  and  the  normalization  condition 
Sfc=l  Dk^lm  T  VO —  1  • 

Theorem  6.3  Let  zq  and  z\  be  1  x  m  and  lxlm  row  vectors,  respectively,  such 
that  the  vector  z  ( zq,z\ )  is  the  unique  positive  solution  to  the  equations 

Then,  the  stationary  probability  vector  y  —  (yo,yi,  ■  ■  ■)  has  the  following  matrix- 
geometric  form 

2/o  =  c_1  20  , 

2/i  =  c-1  z1  ,  (6.4) 

2 Ik  =  2/i  ,  k  =  1,2,-  ••  , 

with 

C  —  Zq  6m  T  Z\  ( Llm  R~)  &lm  •  (6.5) 


In  this  special  case,  the  matrix-geometric  nature  of  the  solution  yields  the 

moments  of  the  queue  length  distribution  in  a  much  simpler  form.  Although  in 

principle,  higher  order  moments  can  also  be  obtained,  only  the  first  two  moments 
are  given  here. 

Corollary  6.4  The  first  two  moments  of  the  stationary  queue  length  distribution  at 
arbitrary  times  are  given  by 

E(Y)  =  yi(Iim  ~  R)~ 2  eim  ,  (6.6a) 

E(Y2)  =  2/1  (Lm  +  R)  {hm  -  R)~3  eim  •  (6.66) 

Proof:  The  following  equalities  can  be  shown  by  direct  computation; 

OO 

£  kRk~x  {Ilm  -  R )2  =  Ilm  (6.7 a) 

k= 1 

OO 

]P  k2 Rk~x  (Jlm  -  Rf  =  Ilm  +  R  (6.76) 

k=  1 
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Equation  (6.6a)  now  follows  from  equations  (6.4),  (6.7a)  and  Lemma  6.1  since 


E(\  )  —  ^  '  ky^cim 

k= l 


OO 

2/i  ^  kRk~1etm 

k=  l 


—  2/l(-Lm  R)  6 Im  • 


Similarly,  the  second  moment  follows  from  equations  (6.4),  (6.7b)  and  Lemma  6.1  . 

n 


Note  for  the  M/M/1  case  that  R  —  p  and  (6.6)  reduces  to  the  well-known 
expressions 


E{Y)  = 


E(Y 2)  = 


P 

l-p' 

P  |  P2 

(l-p)2  +  (l-p)2  • 


NUMERICAL  IMPLEMENTATION 

The  matrix  R  of  equation  (6.2)  can  be  computed  by  using  any  of  the  algorithms 
summarized  in  the  Appendix  by  replacing  the  upper  limit  in  the  summation  for  n 
by  2.  The  comments  following  equation  (5.6)  for  the  computation  of  the  G  matrix 
also  apply  to  these  algorithms.  However,  since  the  matrix  R  is  not  stochastic,  the 
extrapolation  (5.9)  and  the  stopping  criterion  (5.10)  cannot  be  employed  to  the  R 
matrix  calculations  in  general.  Therefore,  in  order  to  obtain  reasonable  accuracy  e 
must  be  chosen  sufficiently  small  and  this  may  require  huge  number  of  iterations, 
especially  for  large  utilizations. 

However,  in  this  quadratic  case,  the  extension  proposed  in  Section  5  can  also  be 
used  in  the  computation  of  the  R  matrix,  in  view  of  the  results  obtained  by  Latouche 
[10].  In  summary,  Latouche’s  paper  establishes  several  useful  relationships  between 
the  matrices  R  and  G  for  the  quadratic  case.  In  particular,  Latouche  shows  that,  if 
the  matrices  R  and  G  are  the  minimal  nonnegative  solutions  to  the  matrix  quadratic 
equations 


R  -  A0  +  RAi  4-  R2A2,  and  G  =  A2  +  AtG  4-  A0G2,  (6.8) 
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then 


R  =  A0(I-A1-  AoG)-1,  and  G  =  (I  -  Ax  -  RA2)~1A2  .  (6.9) 

In  view  of  this  result,  the  matrix  R  can  also  be  computed  first  by  computing 
the  G  matrix  with  any  of  the  algorithms  (A.i)-(A.iv)  of  the  Appendix  together  with 
the  proposed  extrapolation  technique,  and  then  obtaining  R  from  the  first  equation 
in  (6.9).  Note  for  the  MMPP/PH/ 1  queue  that,  the  matrix  R  is  of  dimension  ml. 

An  alternative  approach  is  not  to  use  the  special  structure  of  the  service  time 
distribution  and  solve  an  equation  of  the  form  (1.1)  for  a  G  matrix  using  the  ran¬ 
domization  algorithm  given  by  equation  (2.5).  Although,  the  equation  is  of  order 
N  after  truncation,  the  matrix  G  is  only  of  dimension  m.  However,  the  analysis 
of  the  stationary  probabilities  for  the  MjG / 1  type  queues  is  much  less  transparent 
and  the  results  cannot  be  put  in  as  concise  and  explicit  a  form  as  in  the  matrix- 
geometric  method  [14].  The  numerical  discussion  below  aims  to  shed  some  light  in 
providing  guidelines  in  choosing  the  appropriate  methodology  in  cases  where  both 
are  applicable.  Numerical  results  for  the  M/G/l  approach  corresponds  to  the  RA 
when  used  with  the  stopping  criterion  (5.10),  called  RAea;*r  hereafter. 

The  first-order  computational  requirements  of  the  algorithms  (see  Appendix), 
per  iteration ,  for  the  MMPP/PH /I  model  are 

SS  and  MSS  :  3(/m)3,  MNK  :  8 {Imf,  RA  :  2 Nm3. 

One  would  compare  the  computational  complexities  of  the  matrix-geometric  and  the 
M/G/l  approaches  by  comparing  N  to  l3.  However,  the  results  of  Table  6.1  shows 
that  even  when  l  is  small,  the  M/G/l  approach  yields  far  less  number  of  iterations 
and  gives  much  more  accurate  results  compared  to  any  of  the  algorithms  (A.i)-(A.iv). 
Furthermore,  as  before,  the  number  of  iterations  in  RA  did  not  increase  with  the 
utilization  p  (although  N  slightly  did),  making  the  M/G/l  approach  especially 
superior  at  higher  utilizations.  In  Table  6.1,  the  notation  X  and  Xexir  denote  the 
algorithm  X  when  used  with  the  stopping  criterion  (5.8)  and  (5.10),  respectively. 

The  numbers  in  Table  6.1  for  algorithms  (A.i)-(A.iv)  also  compares  these  algo¬ 
rithms  with  each  other  as  well  as  the  effect  of  the  proposed  extrapolation  method 
on  these  algorithms.  For  further  examples  and  comparisons  the  reader  is  refered  to 
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[5].  It  is  concluded  in  [5]  that  for  small  values  of  N ,  the  MSSe:c<r  is  to  be  prefered 
over  the  MNI(ea:fr  due  to  the  relative  effect  of  the  overhead  in  the  MNK  scheme. 
The  comparison  of  the  first  order  operation  count  per  iteration  given  above  and 
the  number  of  iterations  required  for  convergence  in  each  scheme  in  Table  6.1  also 
supports  this  claim.  Therefore,  the  M/G/l  approach  (RA)  is  compared  only  to 
MSSextr  in  the  following  discussion;  Note  that  for  N  ss  |/3,  the  RA extr  and  the 
MSSextr  have  approximately  the  same  operation  count  per  iteration.  In  Examples 
6. 1.2-6. 1.4,  even  when  /  =  2,  MSSea;tr  and  RAea;tr  required  about  the  same  CPU 
time  per  iteration.  However,  the  number  of  iterations  in  RAea;fr  were  far  less  than 
the  number  of  iterations  in  MSS^,..  RAea;*r  gave  the  same  accuracy  in  the  average 
queue  length,  E(QL ),  with  savings  by  a  factor  of  approximately  25  (resp.  5)  in 
Examples  6. 1.1-6. 1.3  (resp.  in  Example  6.1.4).  In  Example  6.1.5,  although  RAexir 
required  about  3.5  times  more  CPU  time  per  iteration  than  MSSea,<r,  it  resulted 
in  savings  in  the  number  of  iterations  by  a  factor  of  20,  thus  yielding  savings  by  a 
factor  of  approximately  6  in  CPU  time.  In  Examples  6.1.7  and  6.1.8  since  l  =  1, 
the  MSSezfr  required  much  less  CPU  time  (approximately  by  a  factor  of  20),  but 
even  this  is  offset  by  the  huge  difference  in  the  number  of  iterations;  a  factor  of  20 
and  140  in  Example  6.1.7  and  Example  6.1.8,  respectively.  On  the  other  hand,  in 
Example  6.1.6  where  /  =  4,  the  MSSea,tr  required  about  3  times  more  CPU  time  per 
iteration.  Note  also  in  this  example  that  one  gets  the  same  accuracy  in  RAe;r<r  in 
about  20  times  fewer  number  of  iterations,  i.  e.,  savings  by  a  factor  of  60  in  CPU 
time.  Potential  savings  that  can  be  obtained  by  choosing  RAea;ir  over  MSSea;tr  for 
larger  values  of  /  is  no  further  elaborated  but  left  to  the  readers  imagination! 
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TABLE  6.1 


(i)  SSea?^r,  (ii  )  MSS,  (iii)  MSSea;<r,  (iv)  MNK,  (v)  MNKertn  (vi)  RAea;tr. 


#  of  iter.  E(QL) 


o 

bX) 

O 

7 

(i) 

(ii) 

(iii) 

(iv) 

(v) 

(vi) 

(i) 

(ii) 

(iii) 

(iv) 

(v) 

(vi) 

Example  6.1.1. 

m—2,  1-2. 

3 

49 

86 

15 

77 

13 

10 

40.096 

22.972 

29.028 

23.583 

29.162 

26.797 

N=18.  /)=0.768. 

5 

161 

315 

120 

263 

112 

17 

26.905 

26.746 

26.842 

26.756 

26.835 

26.797 

6 

423 

437 

229 

362 

203 

20 

26.815 

26.792 

26.803 

26.793 

26.802 

26.797 

7 

853 

559 

349 

460 

300 

20 

26.799 

26.797 

26.798 

26.797 

26.798 

26.797 

Example  6.2.2. 

m=2.  1=2. 

5 

163 

420 

77 

378 

74 

26 

93.458 

80.952 

91.723 

82.060 

91.585 

89.395 

N=13.  p=0.987. 

7 

655 

1220 

371 

1059 

349 

40 

89.582 

89.293 

89.475 

89.308 

89.466 

89.395 

8 

1431 

1636 

709 

1413 

649 

47 

89.430 

89.385 

89.405 

89.386 

89.403 

89.395 

Example  6.1.3. 

m=2.  1=2. 

5 

100 

85 

50 

51 

47 

9 

6.5371 

6.4897 

6.5099 

6.4915 

6.5085 

6.5000 

N=13.  p=0.875. 

7 

395 

149 

113 

131 

102 

11 

6.5006 

6.4999 

6.5001 

6.4999 

6.5001 

6.5000 

9 

789 

213 

177 

186 

157 

14 

6.5000 

6.5000 

6.5000 

6.5000 

6.5000 

6.5000 

Example  6.1.4. 

m=2.  1=2. 

3 

50 

24 

14 

22 

13 

9 

3.2676 

3.1497 

3.2611 

3.1584 

3.2620 

3.2126 

N=ll.  p= 0.767. 

5 

132 

57 

44 

51 

40 

18 

3.2122 

3.2110 

3.2122 

3.2111 

3.2122 

3.2116 

7 

321 

91 

77 

80 

69 

28 

3.2117 

3.2116 

3.2116 

3.2116 

3.2117 

3.2116 

Example  6.1.5. 

m=2.  1=2. 

5 

114 

93 

54 

83 

51 

9 

7.2496 

7.1974 

7.2187 

7.1988 

7.2172 

7.2083 

O 

CO 

CO 

o 

II 

II 

2 

7 

466 

163 

122 

144 

111 

15 

7.2091 

7.2082 

7.2084 

7.2082 

7.2084 

7.2083 

9 

961 

232 

192 

204 

171 

22 

7.2083 

7.2083 

7.2083 

7.2083 

7.2083 

7.2083 

Example  6.1.6. 

Jl 

<N 

II 

E 

5 

110 

93 

53 

83 

50 

9 

6.9253 

6.8748 

6.8963 

6.8764 

6.8948 

6.8856 

N=35,  p= 0.888. 

7 

440 

163 

122 

143 

110 

16 

6.8855 

6.8855 

6.8857 

6.8855 

6.8857 

6.8856 

9 

668 

232 

192 

173 

140 

22 

6.8857 

6.8856 

6.8856 

6.8856 

6.8856 

6.8856 

Example  6.1.7. 

m=2.  1=1. 

5 

73 

71 

47 

53 

39 

5 

5.0092 

4.9962 

5.0038 

4.9972 

5.0027 

5.0000 

N=28.  p=0.833. 

7 

186 

119 

95 

88 

73 

6 

5.0001 

5.0000 

5.0000 

5.0000 

5.0000 

5.0000 

Example  6.1.8. 

m=2,  1=1. 

5 

78 

259 

59 

207 

54 

5 

51.742 

46.433 

51.224 

47.354 

51.079 

50.000 

N=32.  p=0.980. 

7 

381 

707 

271 

530 

234 

6 

50.068 

49.959 

50.036 

49.991 

50.028 

50.000 

9 

1313 

1169 

707 

860 

554 

8 

50.001 

50.000 

50.004 

50.000 

50.000 

50.000 

-25- 


The  above  discussion  of  the  results  in  Table  6.1  indicates  that  even  when  the 
number  of  service  phases  is  small  there  is  computational  advantage  in  ignoring  the 
special  structure  of  the  service  time  distribution  and  using  the  M/G/ 1  approach. 
It  should  be  reminded  however  that,  as  already  indicated  above,  both  the  waiting 
time  and  the  queue  length  distribution  calculations  are  much  simpler  in  the  matrix- 
geometric  case  once  the  R  matrix  is  obtained.  In  the  M/G/ 1  approach,  as  given 
here  in  Section  3  for  the  MMPP/G/ 1  queue,  these  distributions  are  typically  given 
in  terms  of  the  transform  equations,  whence  Laplace  and  Z  transform  inversions 
are  needed.  Since  the  transform  inversion  methods  require  considerable  CPU  times, 
especially  when  accuracy  is  of  primary  concern,  the  decisions  must  be  based  on  the 
number  of  distribution  points  and  the  degree  of  accuracy  desired  in  these  distribu¬ 
tions.  However,  for  the  MMPP/G/ 1  queue,  Lemmas  3.1  and  4.1  gives  the  moments 
of  the  queue  length  and  the  waiting  time  distributions  recursively  in  a  fraction  of 
the  time  required  to  compute  the  G  matrix  without  having  to  compute  the  inverse 
transforms  and  the  matrices  {An}o°.  Also,  an  alternative  approach  in  which  an  R 
matrix  that  satisfies  equation  (6.4)  but  not  (6.2)  is  obtained  explicitly  in  terms  of 
the  system  parameters  is  presented  in  [4]  for  the  MMPP/ PH /1/K  queue. 
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APPENDIX 

REVIEW  OF  THE  ITERATIVE  ALGORITHMS 


In  this  Appendix,  iterative  algorithms  available  for  finding  the  matrices  G  and 
R  in  equation  (1.1)  are  briefly  surveyed.  Although  the  discussion  is  given  here  for 
the  G  matrix,  these  algorithms  apply  mutatis  mutandis  to  the  calculation  of  the  R 
matrix.  It  is  assumed  that  p  <  1,  or  equivalently,  that  the  matrix  G  is  stochastic. 

(A.i.)  Successive  Substitutions  (SS):  It  is  shown  in  [12]  that  the  sequence  of 
matrices  {G/Jq0  obtained  by 


OO 

Go  =  O  ,  Gk+i  =  Y  AnGnk  ,  k  >  0, 

71=0 


(A.I) 


is  non- decreasing  and  converges  to  the  G  matrix.  The  convergence  rate  of  this  direct 
iterative  scheme  is  R-linear ,  and  can  be  extremely  slow  when  p  is  close  to  one. 

(A.ii.)  Modified  Successive  Substitutions  (MSS):  In  cases  where  the  matrix 
(/  —  A\)  is  well  conditioned,  there  is  computational  advantage  in  using  the  iteration 


G  o  =  0, 


Gk+1  =  (/—  Ax)”1 


A0  +  Y  AnGl 


71=2 


k  >  0.  (A. 2) 


Note  that,  since  the  matrix  A\  is  substochastic,  (/  —  Ai)-1  exists.  It  can  easily  be 
shown  that  the  sequence  of  matrices  {Gfc)o°  obtained  from  MSS  is  componentwise 
larger  then  the  corresponding  matrices  obtained  from  SS.  Therefore,  although  MSS 
requires  a  matrix  inversion  it  converges  in  fewer  number  of  iterations.  It  is  reported 
in  [19]  that  MSS  results  in  about  20%  savings  compared  to  SS. 

(A.iii.)  Newton-Kantorovich  Method  (NK):  The  Gateaux  derivative  [16]  of 
the  mapping  F(X)  =  X  —  AnXn  at  X  is  given  by  the  linear  map 


OO  71  —  1 

[F'(X)]  :U~U-YY  AnXlUXn~1-1 

71=1  1-0 


It  is  shown  in  [13]  that  F'(X)  is  a  non-negative  matrix.  Furthermore,  for  O  <  X  <  G 
(componentwise),  F'(X)  is  an  isotone  so  that  the  operator  F  is  order-convex.  It 
then  follows  by  the  Monotone  Convergence  Theorem  [16,  p.  45]  that  the  Newton- 
Kantorovich  scheme  for  (1.1)  given  by 

Go  =  O  ,  Gfc+1  =Gk-  [F’iGk)}-1  F(Gk),  k  >  0, 


converges  monotonically  to  G.  The  iterates  {Gfc)o°  can  equivalently  be  obtained  by 
the  iteration  [13,  19] 


Go  =  O  ,  Gk+i  =Gk  +  Yk ,  k  >  0,  (A. 3. a) 
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where  Yk  is  the  unique  solution  to  the  linear  system 

oo  n— 1 

Yk  =  -F(Gk)  +  EE'1  nGlkYkGnk-1-1  k>0.  (A.3.b) 

n= 0  1=0 

The  numerical  results  in  Neuts  [13]  indicate  that  despite  the  faster  rate  of  conver¬ 
gence,  the  NK  scheme  is  more  time  consuming  than  the  first  two  schemes,  as  the 
system  in  (A.3.b)  needs  to  be  solved  at  every  iteration.  Therefore,  this  scheme  is 
not  considered  in  this  paper. 

(A.iv.)  Modified  Newton-Kantorovich  Scheme  (MNK):  The  following  mod¬ 
ification  is  proposed  by  Ramaswami  [19]  for  the  NK  scheme  in  order  to  avoid  the 
solution  of  the  linear  system  in  (A.3.b). 

Go  =  O  ,  Gk+ 1  =  Gk  +  Yk,  (AA.a) 

where 

Yk  =  —F(Gk)  +  A\Zk  +  A2{ZkGk  +  GkZk) ,  (AA.b) 

with 

Zk  =  -(I-A1)~1F(Gk)  k>  0.  (A.4.c) 

The  algorithm  (A. 4)  corresponds  to  truncating  the  summation  in  (A.3.b)  at  n  =  2 
and  using  the  estimate  Zk  of  Yk  instead  of  solving  the  linear  system.  It  is  shown  in 
[19]  that  the  MNK  scheme  also  converges  monotonically  to  G. 
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